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ABSTRACT 


This final technical report describes numerical procedures formulated 
for the analysis of three-dimensional (3D) jet mixing problems, as 
incorporated in the computer model, SCIP3D. The overall methodology closely 
parallels that developed in the earlier 2D/axisymmetric jet mixing model, 
SCIPVIS. SCIP3D integrates the 3D parabolized Navier-Stokes (PNS) jet 
mixing equations, cast in mapped cartesian or cylindrical coordinates, 
employing the explicit MacCormack Algorithm. A pressure-split variant of 
this algorithm is employed in subsonic regions with a sublayer approximation 
utilized for treating the streamwise pressure component. SCIP3D contains 
both the ke and kW turbulence models, and employs a two-component mixture 
approach to treat jet exhausts of arbitrary composition. Specialized grid 
procedures are used to adjust the grid growth in accordance with the growth 
of the jet, including a hybrid cartesian/cylindrical grid procedure for 
rectangular jets which moves the hybrid coordinate origin towards the flow 
origin as the jet transitions from a rectangular to circular shape. 
Numerous calculations are presented for rectangular mixing problems, as well 
as for a variety of basic unit problems exhibiting overall capabilities of 
SCIP3D. 



1. INTRODUCTION 


1.1 Program Overview 

This final technical report describes computational methodology 
developed to analyze 3D jet mixing problems using spatial matching 
parabolized Navier-Stokes (PNS) methodology. The program was initiated in 
1981 with preliminary work focused on developing methodology for analyzing 
2D supersonic free jet interaction problems. The SCIPVIS shock-capturing 
model, initially developed to analyze such supersonic problems, was then 
extended to analyze supersonic jets exhausting into subsonic external 
streams. The SCIPVIS methodology was documented in an interim technical 
report, 1 and a number of journal articles were published which describe its 
application to analyzing the detailed multiple-cell shock structure in 
turbulent jets (see refs. 2-5). 

The SCIPVIS model employed explicit numerics and thus could not readily 
analyze wall-bounded jet flows which resolve the details of the near-wall 
boundary layer. To address this problem area, the 2D SPLITP model was 
developed which employs implicit, pressure-split numerics. A number of 
earlier papers and journal articles describe various aspects of the numerics 
in SPLITP (see refs. 6-9). The extension of SCIPVIS and SPLITP based 
numerics to analyze 3D jet mixing problems is the subject of this final 
technical report. 

1.2 Utilization Of SCIPVIS and SPLITP Models In Other Programs 

The SCIPVIS and SPLITP 2D research models developed under this NASA 
program have been extended to analyze a broad spectrum- of problem areas 
under other government sponsored programs. Under Army (MICOM) and Air Force 
(AFWAL) support, SCIPVIS and SPLITP have been extended to analyze 
multi-component/multi-phase flows with nonequilibrium chemistry and 
gas/particle interactions. They serve as components of the latest JANNAF 
Standard Plume Flowfield Model (SPF/3) as schematized in Figure 1, and 
described in refs. 10 and 11, and in several papers presented at the JANNAF 
15th and 16th Plume Technology Meetings. Analogous pairing of SCIPVIS and 
SPLITP PNS numerics has been utilized to analyze hypersonic wake flowfields 
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FIGURE 1. Utilization Of SCIPVIS and SPLITP Models For Zonal PNS Analysis 
Of Rocket Plume Flowfields In JANNAF SPF/3 Standard Plume 
Flowfield Model. 


2 


as described in refs. 11-13. In both the plume and wake models above, 
SCIPVIS and SPLITP are also utilized to analyze the base region employing 
Chapman-Korst type assumptions to achieve closure (see refs. 10 and 14 for 
details) . A wall jet version of SPLITP is a major component of the 
circulation-control airfoil code 1 *- 17 developed jointly by SAIC and AMI 
under NASA /Ames support. 

The most recent extensions of SPLITP and SCIPVIS numerics were made 
under DARPA/Air Force Support to analyze scramjet propulsive flowfields in 
support of the National Aerospace Plane (NASP) program. The NASP version of 
SPLITP, entitled SCORCH, is utilized to analyze scramjet supersonic 
combustor flowfields with tangential or moderately inclined fuel injection 
(refs. 18 and 19). The NASP version of SCIPVIS, entitled SCHNOZ, is 
utilized to analyze scramjet nozzle flowfields including the plume-like 
interaction of the nozzle exhaust flow with the vehicle aerodynamic flow at 
supersonic flight conditions (refs. 20 and 21). These codes have enhanced 
numerical capabilities to deal with strongly combusting flows in comparison 
with the earlier plume/wake versions. SCORCH and SCHNOZ comprise components 
of a complete design-oriented scramjet propulsive flowfield model 21 , 2 * which 
has seen widespread usage by the NASP community. Current efforts are 
focused on unifying SCORCH/SCHNOZ numerics for the multi-zone analysis of 
generalized 2D jet mixing problems with chemistry, utilizing SCORCH implicit 
numerics in near-wall regions, and SCHNOZ explicit numerics in central 
regions . 

1.3 3D Jet Mixing Code, SCIP3D 

The SCIP3D code represents extensions of the 2D methodology in the 
SCIPVIS and SPLITP models to analyze 3D jet mixing problems. SCIP3D is 
based upon the numerics in the SCIPVIS code and employs the explicit 
MacCormack algorithm to spatially integrate the 3D PNS jet mixing equations 
in supersonic flow regions utilizing mapped rectangular or cylindrical 
coordinates. Both the ke and kW two-equation turbulence models are 
incorporated into SCIP3D. For subsonic jet mixing problems, a pressure- 
split approach is employed whereby the streamwise pressure gradient is 
imposed and the cross-flow velocity field is determined via a coupled 
noniterative approach which parallels that employed in SPLITP for highly 
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curved jets.* For mixed problems (supersonic jet into subsonic stream), it 
has been found necessary to adapt the grid to the shape of jet sonic line to 
deal with supersonic/ subsonic coupling methodology in an 'organized' manner. 
For rectangular jets with large aspect ratios, this has required the 
utilization of hybrid rectangular/ cylindrical mapped grids with a moving 
cylindrical origin which ultimately coincides with the jet axis far 
downstream when the jet takes on a circular shape. Numerous test cases are 
described in this report which were performed to check-out the varied run 
options included in SCIP3D. Various aspects of SCIP3D methodology have been 
given in earlier papers and journal articles (see refs. 24-27). This report 
will attempt to synthesize the work already documented, expand upon it where 
necessary, and provide details of recent work not yet reported upon. 
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2. GOVERNING EQUATIONS 


2.1 The 3D PNS Equations 


The conservative form of the 3D PNS 'straight-back' jet mixing equa- 
tions in cartesian (J=0) or cylindrical (J=l) coordinates is listed below: 


9E + 3_ 3F_ 
0x x p J dx 2 



+ H *= V p 


( 1 ) 


In eq. (1), x 4 corresponds to the axial (marching) direction, x 2 is the 
cross-flow direction, and x, is the vertical or radial direction. The 
coordinate system nomenclature for the two systems is summarized in Figure 
2. The vector arrays comprising S’, F, S’ and ff are as follows: 


pU 1 


pV 


pW 

aP+pU* 


pUV 


pUW 

pUV 

F = 

P+pV 2 

G = 

pVW 

pUW 


pVW 


P+pW 2 

pUH 


pVH 


pWH 

pUd 


pvZ 


pV6 

, __ 






— 1 

JpW 

r^d-oiap/ax^jpuw 


0 

v u 

H = — 

2JpVW 


V v 

r J 

J(W 2 -V 2 ) 

JpWH 

JpW? 


1 

1 


The vector represents the parabolized stress and transport terms. In 

the above equations, U, V, and W are the axial, crossflow and radial 
velocity components, p is the density, P is the pressure, H is the total 
enthalpy, and Z is the species parameter to be described below. The 
parameter, a, is a pressure-split parameter used to distinguish between 
supersonic (a=l) and subsonic (a=0) run modes. 

2.2 Thermodynamics 


The jet mixing problems considered assume that the jet and external 
streams are each of uniform composition. For nonreacting (chemically 
frozen) situations, the species parameter, Z, describes the local mixture 
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composition; viz.. 


z = 


a. 

l 



a. 

1 


E 



( 2 ) 


where is the mass fraction of the ith species and J and E represent the 
constant values of in the unmixed jet and external streams. The static 
enthalpy is given by: 


h(?,T) = Ch J (T) - h E (T)}£ + h E (T) 


(3) 


where : 


h J,E “ J {a i h i (T)) J,l 


Then the specific heat ratio, y(Z, T) is given by: 


C (d,T)W(d)/R 


C p (<5,T)W(«5) /R 0 -l 


(4) 


where the specific heat, C p , is given by: 


V*- T > - Jr 


(5) 


and the molecular weight, W, is given by: 


W(*» = [ ( -1 ) * + J_ ] 

J E E 


—i 


( 6 ) 


2.3 Parabolized Stress and Transport Terms 

Eliminating all terms containing streamwise derivatives and using the 
Boussinesq-type approximation, the parabolized stress terms in the cartesian 
and cylindrical systems are written as follows: 
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Cartesian Coordinates ^ = x, x 2 » y, x g = z) 


v - 9 t 9V \ j. 9 t 9V \ 

V U 3y (fl 3y ) 3z (fl 3z ) 


„ 3 ,4 3 V* , 3 t 3V. 3 / 2 3W. . 3 - 3W. 

V V “ 3y ( 3 ** 3y ) dz** 1 3z ) 3y ( 3 * 3z* 3z (fl 3y ) 


„ 3 /4 3W. . 3 , 3Wv 3 / 2 3V. . 3 , 3V. 

V W 3z ( 3 11 3z ) 3y (t * 3y ) 3z 3 M 3y } 3y i]l 3z 


H 


L. (? “) + 1. (? “) 

3y M 3y' 3z M 3z' 


3 / ~ .3 / ~ 3tf% 

Sy ^ 5y> + ^ 3J> 


Cylindrical Coordinates (x^^ = x, x^ = x^ = r) 


„ i 3 , i 30. . a , ao. 

v u = F 3? ^ F dF } sF ^ sF ) 


u a_ 

3 r 3r 


v _ 1 <L 

H r dr 


1_ a_ [ r «„ (9V + i 4-1)1 + i ® L <i I 2} 

r* 3r [ r 11 l 3r r 3F r ; J r dF l 3 r dF 

[rW] ) + 2|i £ j 

li^L. 

3 H 3r 

T ,3V , 1 3W . V 1 1 \ 4 1 3V 

1/ ( dF F dF )_ ^ fJ " F “[? F aF " 


(rW) " 1 11 If _ 2vl w 1 + 


■] 


1 3 
r dF 


J I ( rW ) + 2W/r 
3 r dr 


] 


(7a) 

(7b) 

(7c) 

(7d) 

(7e) 

(8a) 

(8b) 

(8c) 
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(8d) 



d 

n 


r 1 3H \ 

<|l r 


a aH x 

3r ‘o aF> 




3 3<5. 

1? «7> 


(8e) 


In the above stress/transport term expressions, p represents the 'effective 
viscosity' (= p.^ + p^) and p is the effective transport parameter for scalar 
diffusion (« p^yPr + p fc /Pr fc ) based on the assumption that heat and mass 
diffuse at the same rate (viz., Lewis number of unity assumed). 

2.4 Two-Equation Turbulence Models 

Both the ke 1 * and kW*» two-equation turbulence models are utilized in 
SCIP3D. The standard coefficients and constants are employed as described 
in the earlier SCIPVIS report,. 1 Using tensor notation (repeated subscripts 
indicating summation), the ke model solves the following equations for k and 



(pVjk) = 



H 

a k 



+ P - 8 


(9a) 


fe: ‘'V ■ fe: < r m ht > + c i f 

J J 8 j 


where : 


dU, 3U. 3U 


wu • ww . 

p = u ( — i + — =■ ) 

“ “t < 3x ± 3Xj 


IT" h m 


(9b) 


( 10 ) 


The dominant cartesian parabolized terms in the vorticity, <o, for flows 
where U >> V, W are given by: 


i 


w 


"t 


<§£>* + <£>* + 0 



(id 
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The dominant cylindrical parabolized terms in the vorticity, o>, are given 
by: 


<» = p. I 


, i au 

( 2 U 


) 1 + 


( 

1 dr 


) + 0 C [ 


atv.w] 
m~FT 1 


The turbulent viscosity, p^., is defined by: 


p fc = C^pk /e 


( 12 ) 


The kW model solves the following equations for k and W: 


(p V° 


® ( Hi ) + p - P c kw 1/2 

3x j °k 9x j " D 


(13a) 


k; (pV j w) 


3 

8x 


J ' 5 KJ > + T= - V" 3/2 * Vt""' 


(13b) 


The dominant cartesian parabolized terms in the vorticity expression, (V«d) , 
are enumerated below: 


m i [\au a*u/ . , au au a*u 
( V „), = |^<3y JyT> 


a*u + ,au a»u } 

dydz dz 3y3z ; 


f au a*u > + 2 22 a * u 

l 3y az3y ; 3y dz 3z a 3z3y 


+ (22. £!£)* 1 

az az* J 


(14) 
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The dominant cylindrical parabolized terms in the vorticity expression, 
2 

(Vw) , are given by: 


(v.)* = + 2 to •!£ + <1 SI) L. (I W)1 

r* r |_ v 3r 3r* > 3d' 3r > W j 

+ 1 _ r ( 9 £v a^u 1 au a_ ,1 aikl* 

w a [_ 3r 3r* r 3? 3r r 3© J 

+ iTlE I 3 A . i au i a ,1 au/l* 

u 1 [ar r d? 3r r 3< r 3? r 3© J 


(15) 


The turbulent viscosity, n fc , is defined by: 

H t = pk/W 1/2 (16) 

2.5 Happed Vectorized Equations 


The mean flow and turbulence model equations are solved in a mapped 
computational domain defined by the geometric transformation: 


5 = x i 


11 X 2 /x 2,U 

4 = (X 3- X 3,L )/(X 3,IT X 3,L ) 


(17) 


In this transformation, L and U designate bounding surfaces of the 
computation. In the present application of this mapping for rectangular 
problems, the bounding computational surfaces are assumed to have no cross 
curvature i.e. y^ = y^fx), \ \j = \ ^(x)* With this assumption, the 

cartesian equations can be cast in the following form: 


3E 3F 
a$ an 



(18) 
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where : 


E = E 

F “ b / " a y E 
G - b 2 G - a 2 E 
H = H + ( ta y ]^ + [a E l )E 

and the transformation parameters a , a , b and b are defined by: 

y z y z 

a y = h(y u ) x /y u 

a 2 = l(l-«)[^] x + t[z 0 ix 1/(z u - 

V 1/y u 

V 1/(z u “ Z L ) 

For free jet problems, y^ and z^ y must fully encompass the viscous 
jet. This is accomplished using growth rule formulations of the form: 


dz u _ ( JMAX-1) .W. Z u ( af. 

dx - (J*-l) fj «Zj MAX 


(19) 


which parallel those utilized in SCIPVIS 1 for 2D underexpanded jets, applied 
on the two symmetry planes (y=0 and Z=0) as exhibited in Figure 3. In 
Figure 3, J* corresponds to the viscous dividing streamline position on the 
symmetry plane, y * 0, ascertained by monitoring the 4 profile at each step 
{4 - 4 * , where 4* = .5 until the mixing reaches the axis); JMAX is the index 
of the upper boundary point; W/U is the streamline inclination at J*; and, 
df/dz is the outer edge gradient of the parameter f (taken to be the maximum 
of the streamwise velocity and species parameter) . This boundary growth 
formulation provides for exact alignment of the viscous dividing streamline 
position with the mapped coordinate line, c = constant in the initial 
expansion region of the jet. Note that the viscous part of the growth 
relation (the df/dz term) is not activated until the jet shear layer reaches 
the upper boundary, and, that boundary growth is enforced to always be 
positive. Application of this boundary growth formulation on both symmetry 
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planes to yield the variation of Zy(x) and y 0 (x) yields a mapped grid 
network which initially captures the 3D jet induced bow shock, and encloses 
the viscous jet in a manner which automatically makes the transition from a 
rectangular to a square type of grid in accordance with the respective 
vertical and transverse gradients in the flow. 

For cylindrical problems, the cross-flow boundaries are always constant 
(i.e. x 2 y = 4y = constant) while the radial boundaries can have both axial 
and cross-curvature (i.e. *3 j = ^ [j = ^ j (x,4). With these stipu- 
lations, the cylindrical equations can be cast in the same form as eq. (20) 
where now: 


E = E (20) 

F = b,P 

G = rb r G " a r,«S ¥ 

H = H + (a ,)$F - G 
r, c 

- <V x >E/(b r r> - Vx r(?/r,t 

and the transformation parameters a ,, a . b. and b„ are defined by: 

r,0 r,x v r 


a r ,6 = r(1 “« )r U + « r U <5 1/(r U- r L ) 
a r,x * t(1 "« )r L + t r „ 1/(r u- r L ) 

X A 

b , - 

b r ' 1/(p lT r L ) 


For rectangular jet problems with a large aspect ratio (> 2/1), a 
hybrid mapped rectangular/cylindrical grid (Figure 4) would be utilized. As 
the solution evolves and the jet becomes more circular, the cylindrical 
origin moves to the flow origin. The procedures entailed in this moving 
origin methodology will be described in a subsequent section of this report. 
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FIGOBE 4. Hybrid Rectangular/Circular Mapping In Cross-Flow Plane For 
Generalized Rectangular Jet Problems. 


5 - 




3. NUMERICAL PROCEDURES FOR AHALYZIHG SUPERSONIC 
MIXTHG PROBLEMS ABB NUMERICAL STUDIES 


3.1 Interior Point Algorithm 

In present versions of SCIP3D, the grid points are equally spaced in 

both the 11 and £ directions. All convective terms are represented using 

one-sided (two point) difference expressions; diffusive terms are 

represented using central difference expressions. The solution at interior 

grid points is obtained using a spatial marching, rotating version of the 

MacCormack algorithm (see ref. 30) . For the grid point, I, J (see Figure 

3 ) , this two-step algorithm takes the following form in advancing the 

solution from 5 to ? + A£ (~ denotes predictor level values; - denotes 

corrector level values at £ + A §; and, e and are alternated between 0 

y z 

and 1 at odd/even steps to provide a nonpreferential (rotating) treatment of 
wave/ convective processes) . 


Predictor Step 


E(I.J) - E(I.J) - |£ [d-« y )F I+1>J - <1 - VW 

" At I(1 ~'z )G I,J+l ■ (1_2 'z )G I,J “ 'z G I.J-l J ■ H I,J A5 
+ V F (I,J)A5 


(20a) 


Corrector Step 

e(i 


• J ) - I ( E<I - J > * + u ' 2 , y )F I,J 

Uy-DFj.j jl - Ie z G I<J+1 + ,J + (e z _1) °I, J-l 1 

- h^jAS + v p (i,j)A5 | 


(20b) 
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where, for f = D, V p takes the form: 


b 2 ^ 

V p (I,J) = a ^ { t 1 I I f ( I+1 * J > - - f(I-l.J)] } 


( 21 ) 


b 2 


+ { 4j[f(I»J+l) ~ f(I»J)3 - Hj[f(I,J) - f(I,J-l)J } 


After both the predictor and corrector steps, the conservation array, 
E(I,J), is decoded using standard supersonic flow decode procedures (see 
ref. 1). 

The algorithm utilized is fully explicit and the step-size is limited 
by both hyperbolic (CFL) and parabolic stability constraints. These are 
combined in a 'parallel-resistor' manner (viz., A£ = (l/A5 Hyp + l/A£p AR )- 2 ) 
with the CFL constraint satisfied using nonlinear reference plane 
characteristic intersections. For the fully supersonic jet mixing problems 
under consideration, the hyperbolic CFL constraint effectively controls the 
allowable step-size, with the parabolic criterion having little impact on 
the allowable step-size. 

3.2 Grid Distribution and Step-Size Control 

Consider the rectangular jet mixing problem with mild (e.g., 2/1) 
aspect ratio depicted in Figure 5. A simplistic hybrid grid (drawn with 
extremely crude grid definition) could be utilized to initiate the 
calculation. As the jet contour takes on an elliptical appearance, the 
origin would be moved to the left (e.g., rectangular planes are eliminated) 
until far downstream (when the cross-section is circular), the cylindrical 
origin coincides with the flow axis of symmetry. Examples illustrating this 
'moving origin' methodology will be exhibited. 
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FIGURE 5. Initial Hybrid Grid For Analysis Of Rectangular Jet Mixing 
Problem. 
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The coupling of cartesian and cylindrical grids is relatively 
straightforward, but issues must be addressed in treating grid points near 
the cylindrical grid origin. Presently, a 'brute-force' approach is 
employed which utilizes a number of overlap cartesian planes (Figure 6), 
with properties interpolated from the cylindrical grid solution. Using 
these overlap planes, the rectangular portion of the flowfield is integrated 
first at the predictor level to establish boundary conditions for the 

cylindrical region, which is then integrated at the predictor level. The 
rectangular portion is then integrated for the corrector step, followed by 
the integration of the cylindrical region for the corrector step. Only one 
solution (rectangular or cylindrical) need be kept in core storage at a 
time. (Note that for very high aspect ratio problems, core storage can be 
minimized by utilizing several blocks of rectangular grids, each solved 

independently utilizing overlap planes for coupling) . Note that at the 

hybrid origin, W (=0) and V in the cartesian system are single valued but 
in cylindrical coordinates, V = 0 and W = V cos6. Specialized procedures 
must be employed to avoid step-size limitations associated with cross-flow 
differentiation in cylindrical coordinates in the vicinity of the origin due 
to the small physical spacing between grid points. Again, the present 

approach is 'brute-force' (e.g., only selective grid points are calculated 
in accordance with the ratio of rAd to r u A(6) . The hybrid grid has only 
recently been incorporated into the code -and further exploratory work 
appears warranted. 

The 'shear layer' initialization for jet interaction problems utilized 
in SCIP3D is directly analogous to that employed in previous overlaid 
viscous/ inviscid jet models.* 1 , ** The initialization entails: 

(a) determining initial mean flow shear layer profiles centered about 
the jet slipstream location a small distance downstream of the 
nozzle exit plane (see refs. 31 and 32 for details); and, 

(b) determining initial turbulence parameter (ke or kW) profiles based 
on an equilibrium turbulence hypothesis and the use of a mixing 
length turbulence model approximation. 
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3.3 Boundary Point Procedures 


The following boundary surface types are analyzed in the supersonic 
version of the SCIP3D model: 

(1) solid surfaces (PNS with slip [no BL] or Euler limit); 

(2) symmetry planes; and, 

(3) jet mixing layer outer edge. 

The treatment of symmetry planes, using appropriate odd/even reflection 
conditions, employs standard techniques. Solid surfaces are treated using a 
variant of the Abbett wave-correction procedure** which incorporates the 
surface entropy approach and corner point methodology of Dash and Del- 
Guidice* 0 . Jet computational boundaries are configured to fall outside the 
viscous jet boundary using the growth rate formulation of equation (19) as 
schematized in Figure 3. For nonuniform supersonic flow surrounding the 
jet, the external flow can be calculated concurrently using the present 
shock-capturing algorithm in the inviscid limit. This is accomplished by 
extending the numerical domain outward sufficiently to capture the bounding 
outer shock using an estimate of the shock growth with a safety factor. For 
uniform surrounding flow, the external flowfield can be approximated using a 
linearized pressure/ flow-deflection relation. In applying this approach, 
the nearfield portion of the jet is solved with the grid encompassing the 
captured jet induced bow shock (Figure 3) as discussed previously. The jet 
boundary is calculated using the Abbett wave-correction procedure as will be 
summarized below. 

The approach taken for jets exhausting into uniform supersonic streams 
assumes that the jet has quarter-plane symmetry and hence, the free stream 
velocity, D, is aligned with the x coordinate direction. Consider the 
calculation for upper jet boundary points on the surface Zg(x) (Figure 7). 
The boundary points, J = JMAX, are first calculated using the interior point 
predictor/corrector algorithm with backward/ backward predictor/corrector 
convective differences in the z direction and standard alternating 
predictor/ corrector differences in the y direction. This yields properties 
0*, V*, P*, M*, etc. which do not satisfy the external flow solution. The 
external flow pressure, Pg, is given by the linearized relation: 
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•22 


P E “ P o 




1 + 


(e*+e*) 


(22) 


where 0* is the flow angle based on the predictor/ corrector solution, viz.: 


0 * - S 


tan- 1 


(V* a -W * a ) 1/2 

0* 


(23) 


(S = 1 if 5* . Tp > 0; S «= -1 if Q* . r p < 0; r p is the position vector to 
the boundary point in the plane x «= const; Q* is the velocity vector, Q* *» 
0*i + V*i y + V*± z ) and 0' is a correction to this angle (obtained via 
perturbing the velocity in the Tp direction only) . The corrected boundary 
pressure is given by an isentropic simple wave rotation from P*, 9* 

conditions, viz.: 


p c - p* 


_ yW 


4 


M a -1 


(24) 


Solving eqs. (22) and (24) simultaneously with Pg set equal to Pg 
yields the corrected boundary pressure and flow angle. The total enthalpy 
remains invariant (H^, = H^,*) and the density is expanded isentropically (pg 
* p*( Pg/P*) 1 ^ 2 ). Hence, the magnitude of the corrected velocity is given by 
(for a perfect gas) : 


2y P r 1/2 

<2 V^’ 


(25) 


and the corrected velocity vector can be written: 
Q c - Q* + q' - D*i x + (V*+v')i y + (W*+w')i z 


(26) 
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(27) 


where (since q* * q'i p ) : 

P 


v ' = q ' 


i 

r p 


z 

w' = q' ( -2- ) 

r p 


Substituting the expressions for v' and w* in eq. (27) into eq. (26), and 
setting the magnitude of | | to the value given by eq, 25 yields a 

quadratic expression for the magnitude of q’; the sign corresponds to that 

of e' (i.e., q - | q' k if 6» > 0; q' - - | q' k if 0' < 0). 

P P 

In applying the Abbett wave-correction procedure at solid surfaces, the 
same predictor/corrector sequence discussed above is utilized. Thus, for an 
upper wall point, z u = F(x,y), the MacCormack algorithm would be employed 
with backward/backward differences in the z direction and alternating 
differences in the y. The resultant velocity vector will not satisfy the 
boundary condition V. 'n' = 0 and must be rotated through the angle 0' given 
by Q* . 'n s « | Q* | sin&* yielding 

-U*F -V*F +W* 

0' = sin- 1 { 2L_ — J (28) 

(1+F *+F *) ' Q* 
x y 


with the corrected pressure, P , obtained via the relation 

c 


In (P /?*) = +yM*e'/(M*-l) l/ * 


(29) 


The corrected density, p , is obtained from the isentropic relation 

c 


p c /p l - p */p l 


(30) 


and the corrected velocity, Q , is obtained from 

c 


7 1 2 

H* = — _ _£ + i Q 
T" 1 P c 2 c 


(31) 
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Then, the corrected velocity components are given by: 


/\ /\ /\ /\ 

Q c * " °c V'cVVs 


where 


/s 

t 


Q, - <&„ . ^ ) £ 
5 - (5, . n' ) 'n 


( 32 ) 


(33) 


For corner points with no cross-curvature (i.e. Z u = F(x), y u = G(x)), the 
concurrent satisfaction of V . In' = 0 on both surfaces is obtained by apply- 
ing eq. (29) additively for rotations through the angles where: 


tanO^ 



(34a) 


tanO ' = G 
2 



(34b) 


3.4 Corner Problem Test Case 


The wave/ shock-capturing capabilities of SCIP3D were first checked out 
by the performance of fundamental 3D corner flow interaction problems. 
Calculations were performed for: 

(1) a double expansion corner; 

(2) an expansion/ compression corner; and, 

(3) a double compression corner. 

In these test calculations, waves generated by two perpendicular walls 
interact with each other. Reference 34 provides details of these elementary 
cases. These calculations were performed starting with uniform supersonic 
flow, discontinuously turning the lower and side walls at x = 0. A mapped 
rectangular domain was implemented for all the calculations employing 
quarter-plane symmetry. The calculations were terminated before the wave 
systems reached the symmetry planes. 
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Boundary conditions at the walls were analyzed using a variant of the 
Abbett wave correction technique modified to employ the entropy variation 
along the grid line one grid interval away from the wall (see refs. 30, 34, 
35). As discussed in Section 3.3: (1) wall point properties are first 
calculated using the MacCormack algorithm with differences evaluated 'into 
the flow' in both the predictor and corrector steps; and (2) the pressure 
and flow angles are then 'corrected ' via a rotation in the plane containing 
the wall normal and predicted velocity vector so as to satisfy the boundary 
condition V . / n' * 0. At the corner point, the Abbett technique is formu- 
lated whereby the properties at the corner are evaluated based on: (1) 
application of the MacCormack algorithm using both Y and Z derivatives 'into 
the flow' on both predictor and corrector steps; (2) a dual correction 
procedure to satisfy V . 'n' = 0 on both intersecting surfaces yielding the 
total turning angle, A» T , i.e., 

AO^, = A9^ + A&2 


where AOj and A0 2 are the turning corrections required on the two inter- 
secting surfaces; and, (3) evaluation of the pressure via the relation 


ln(P c o/P M ) = + 




■A&T 


l/ Rr ~ rT 

T M 


(35) 


where P„ is the value predicted by the MacCormack algorithm and P is the 

M C 

corrected value . This approach has been found to work extremely well . 
Results for the double compression corner test case are described below. 


In this Mach 3.17 case, both the lower wall and side wall were turned 
up by 12.5° generating two intersecting shock waves, as exhibited in the 
insert of Figure 8. The predicted pressure contours after 18 steps are 
shown in Figure 8 and exhibit the classic behavior obtained in experiments 
and using conically invariant marching procedures. The streamwise corner 
pressure variation is shown in Figure 9 and 'settles down' after several 
steps to the correct asymptotic value observed experimentally. The initial 
oscillations occur due to capturing a strong corner shock (12 1 /* 0 + 12 V* 0 
of abrupt compression) with no initial grid resolution. 
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FIGURE 8. 


Pressure Contours For Double Compression Corner 


Pressure ( P/ 



3.5 Square Converging-Diverging Kozzle Problem 


The next sample calculation was performed for a square converging- 
diverging nozzle (Figure 10) which encompasses the double-expansion and 
double-compression wave interactions discussed in Section 3.4. Complete 
details of this flowfield are given in ref. 34 which describes results 
utilizing the earlier BIGMAC 3D nozzle model of Dash and DelGuidice. The 
initial Mach number for this case was 2.94 and the initial pressure was 0.4 
atm. The starting profile was uniform. Figure 11 exhibits predicted 
temperature contours on the symmetry planes (quarter plane symmetry was 
assumed) and pressure contours at the station x » 15 showing the occurrence 
of dual plane shock interactions. Figure 12 exhibits the predicted stream- 
wise variation of Mach number along the axis and comparisons with the 
results of Chang* 4 (which terminated at the shock). Figure 13 exhibits the 
streamwise variation of pressure along the axis and comparisons with results 
of the GIM code* 7 . Figure 14 exhibits the predicted streamwise pressure 
variation along a wall/axis juncture. The results compare favorable with 
those of the GIM code and with previous BIGMAC shock-capturing results* 4 
(not shown) . The predicted streamwise pressure variation along the corner 
is shown in Figure 15. 

3.6 Balanced Pressure Square Jet Problem 


This calculation was performed for a hot Mach 3 jet (T t = 2000 °K) 

V 

exhausting into a colder Mach 3 uniform external stream (T E = 500 °K) - the 
resultant jet velocity ratio, Uj/Ug, was 2/1. The calculation was performed 
with a 21x21 square grid utilizing the 'uniform jet ' startline procedure 
described in Section 3.1, to provide initial shear layer profiles. Of 
Interest is the transition from square to circular type contours. The 
predicted outer jet boundary contours (6 = 0.01) shown in Figure 16 for 
several x stations 16 represents an inert species, viz., 6 - 1 in the 
unmixed jet and equals 0 in the external streams; for Pr *= 1, i also 
represents the nondimens ional axial velocity distributions, i.e., 6 = 

(tf-Ug)/(Uj-Ug) for constant pressure mixing] clearly exhibit this 
transition. At x/rj ~ 20, the boundary contour still has a distinct 
noncircular appearance while at x/r^ ~ 40, the contour starts to appear 
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FIGURE 10. Square Converging-Diverging Nozzle. 
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IGOftB 11. Temperature Contours On Symmetry Plane for Square CD Nozzle 
Problem. 
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FIGURE 12. Streamwise Variation Of Mach Number Along The Axis For Square 
Nozzle Problem. 
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FIGURE 15. Streamwise Variation Of Pressure Along The Corner For Square 
Nozzle Problem. 
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FIGURE 16. Jet Boundary Contours For Balanced Pressure Square Jet 

™~ “ “““ Problem, ke Turbulence Model Utilized. 
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circular. At x/r^ ~ 80, the outer contour is essentially circular as is the 
entire jet solution (Figure 17). 

We note that at x/rj ~ 20, the peak turbulence levels (Figure 18) occur 
on the two symmetry planes and are well off axis, characteristic of shear 
layer-like behavior, while further downstream, peak levels are centered 
about the jet axis. The predicted variation of the turbulent kinetic 
energy, k, and the species parameter, 4 , along the jet axis are exhibited in 
Figure 19 and are compared with SCIPVIS 1 predictions for the 'corresponding' 
circular jet (the square jet was initially 2x2; the circular jet had an 
initial diameter of 2) . The square jet is predicted to mix somewhat faster 
than the circular jet (which was also predicted using 21 grid points and the 
same turbulence model) . If we perform the comparison for jets of equal mass 
flux on a scaled basis (viz., k, 4 vs. x/nA), the circular jet would mix 
faster than the square jet. Figure 20 compares the predicted axial species 
decay along the jet centerline using the ke and kW turbulence models. The 
ke model is seen to predict substantially faster mixing in accord with past 
observations for axisymmetric jet problems (see ref. 1). Figure 21 depicts 
predicted species ( 4 ) and temperature (T) contours in the plane of symmetry 
(y * 0 or z = 0) obtained using the ke turbulence model. 

3.7 Pnderexpanded Square Jet Problem 


The same square jet problem described above was repeated with the jet 
exhaust pressure boosted to 2 atm (Pj/P E = 2/1). The calculation was per- 
formed in both rectangular and cylindrical coordinates. 

3.7.1 Rectangular Coordinate Solution 

The rectangular coordinate calculation was performed with a 21 x 21 
rectangular grid initially spread evenly over the quarter plane domain 0 < y 
< 2 and 0 < z < 2. The calculation initiates with a Prandtl-Meyer fan 
propagating inward from the nozzle trailing edge lip to the axis, and with a 
plume induced bow shock propagating outward. Figure 22 exhibits pressure 
contours for this case on the symmetry planes (y = 0 or z = 0) with the 
extent of the shear layer (.1 < 4 < .9) superimposed. The symmetry plane 
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FIGURE 17. 


6 Contours At x/r. = 80 For Balanced Pressure Square Jet 
Problem, ke Turbuldnce Model Utilized. 
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FIGURE 18. Turbulent Kinetic Energy Contours At x/r = 20 For Balanced 
Pressure Square Jet Problem, ke Turbulence Model Utilized. 
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FIGURE 19 . Variation Of Turbulent Kinetic Energy (XK) And Species Para- 
meter (PHI) Along Axis Of Symmetry For Square (SCIP3D) And 
Circular (SCIPVIS) Jets - kW Turbulence Model Utilized. 
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FIGURE 20. Comparison Of ke And kW Turbulence Model Results For Species 
Centerline Decay For Square Jet Problem. 
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FIGURE 21. Predicted Species (PHI) and Temperature (T) Contours In Plane 
Of Symmetry For Square Jet Problem - ke Model Results. 
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plot closely resembles the analogous 2D/axisymmetric jet solution plot. 
However, contours in the cross-flow planes, x = constant, clearly exhibit 
the 3D characteristics of this flow. 

Figure 23 exhibits several pressure contours in the cross-flow planes 
at stations 0 < x < 10. The square startline solution (x ~ 0) just down- 
stream of the cowl lip depicts the Inward expansion and outward compression 
process quite distinctly. The pattern at x = 2 shows the complex effects of 
the interaction of the initially planar/perpendicular wave fronts. On the 
symmetry planes y = z = 0, the pattern is 2D like but on the 45° plane (y = 
z), the 3D interactions of the waves are quite evident (e.g. the contour of 
minimum pressure centered about y = z ~ .75 results from the interaction of 
the two side expansions, etc) . The wave fronts of the barrel shock can be 
distinguished at x = 6 and the pattern is highly three-dimensional. The 
reflected wave fronts are distinguishable at x = 8 and x = 10. 

Figure 24 depicts 6 contours at these same axial stations showing the 
evolution of the jet shear layer contour from square to rounded shapes. The 
initially square contours take on a somewhat irregular shape initially in 
response to the 3D convective cross-flow velocities generated by the 
pressure field. At x = 10, the contours again become smooth with rounding 
about the 45° plane, and downstream, will take on rounded shapes as in the 
balanced pressure case. The outer jet contours (i = .1) for the first 
several radii are compared on a single plot in Figure 25 for both the near- 
field (0 < x < 10) solution and further downstream to x = 40. 

The SCIP3D predicted nearfield centerline pressure variation is 
compared with SCIPVIS predictions for 'corresponding' planar (2D) and 
circular (AXI) jets (Figure 2.6). The initial square jet pressure decay 
predicted by SCIP3D is bounded by the planar and circular SCIPVIS pre- 
dictions. The square jet recompression initiates downstream of the circular 
case and the pressure variation is more rapid with an overshoot occurring at 
x ~ 8 that is not evidenced in the circular or planar situations. 
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FIGURE 26 . Comparison Of SCIP3D Predicted Centerline Pressure Variation 
For Underexpanded Square Jet Problem With SCIPVIS Predictions 
For Analogous Planar (2D) and Axisymmetric Jet Problems. 


i 


-48- 



3.7.2 Cylindrical Coordinate Solution 


In the rectangular coordinate solution, the boundaries z u (x) and y u (x) 
moved outward to encompass the viscous jet as per the discussion of Section 
2.5 and the Figure 3 schematic. In cylindrical coordinates, the outer 
boundary, r u (x,<6), adapts to the shape of the jet contour. Figure 27 
depicts the cylindrical grid at x = 0 and x = 10 employed for this problem 
(21 points in the r direction and 12 points in the i direction were 
utilized) . Figure 28 exhibits pressure contours in the symmetry planes (t = 
0 and 7T/2) which should be compared with those of the rectangular solution 
(Figure 28) . A comparison of predicted centerline pressures along the jet 
axis, and, predicted jet outer boundary variations in the symmetry place are 
exhibited in Figure 29 for the cartesian and cylindrical solutions. The jet 
boundary variations are seen to be essentially identical. The axis pressure 
variations are comparable for 0 < x < 8 but have some departures downstream, 
possibly due to the different grid resolution in the cross-flow direction 
which could affect the near-axis solution in the vicinity of shock re- 
flection. 

A comparison of i contours at x/r^ = 10 (Figure 30) for the cartesian 
grid and cylindrical grid solutions shows that the mixing characteristics of 
the two solutions is essentially grid independent. Cross-flow velocity 
vectors for the cylindrical grid solution at x/r^ =2, 4, 6, 8 and 10 are 
exhibited in Figures 31A, B, C, D, and E respectively, and clearly exhibit 
the wave effects. At x/ r^ = 2 and 4 the plume is expanding and the cross- 
flow is small relative to the grid lines, 4 * constant. At x/r^ = 6, the 
flow is recompressing and a barrel shock is propagating towards the axis 
(see also, the pressure contours of Figure 23). The cross-flow is seen to 
be highly three-dimensional in the vicinity of the barrel shock. At x “ 8 
the barrel shock has reflected from the axis, and, at x * 10, the reflected 
shock has propagated back into the flow turning it parallel to the axis. 

3.8 Balanced Pressure 4/1 Aspect Ratio Rectangular Jet Problem 

This quarter plane symmetry calculation was performed for an exhaust 
emanating from a rectangular nozzle (0<y<2; 0<z<l/2) with the 
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FIGURE 27. Adaptive Cylindrical Grid; Grid Distribution At x * 0 And 10 
For Underexpanded Square Jet Problem. 
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initial 21x21 grid encompassing the domain 0 < y < 2.5 and 0 < z < 1. The 
computational boundaries Y u (x) and z u (x) adapted to the growth of the jet as 
described in Section 2.5. The calculation was run out to x = 120 using the 
ke turbulence model. The exhaust and external flow conditions were the same 
as for the previous balanced-pressure square jet case, namely Mj « 3, Tj = 
2000 °K and Mg = 3, Tg = 500 °K. 

The nearfield structure of the rectangular jet flow is exhibited by the 
contours of the velocity, U; the species parameter, 4 ; the temperature, T; 
and, the turbulent kinetic energy, k; in the symmetry planes y * 0 and z = 0 
as exhibited in Figure 32-35. Cross-flow contours of 4 at x ■ 5, 10, 20 and 
40 (Figure 36) show the nearfield shape change from a rectangular to an 
elliptical cross-section. Cross-flow contours of k at ■ 10, 20 and 30 
(Figure 37) show peak values of turbulence to be off-axis and initially 
dominated by the stronger mixing in the vertical (z) direction. 

The overall jet structure (nearfield/farfield) is exhibited in the 4 
symmetry plane contours of Figure 38 for 0 < x < 120, and, the 4 cross-flow 
contours at x = 60, 80 and 100 (Figure 39) . The centerline variations of 4 , 
T, and k for 0 < x < 120 are exhibited in Figures 40 - 42. 

3.9 Underexpanded 4/1 Aspect Ratio Rectangular Jet Problem 

The above problem was repeated with the jet pressure boosted to 2 atm, 
and all other conditions remaining the same. The pressure contours in the 
two symmetry planes (y * 0 and z * 0) are exhibited in Figure 43 with the 
shear layer (.1 < 4 < .0) superimposed. Jet contours ( 4 = .1) at x ■ 0, 10, 
20 and 40 (Figure 44) exhibit a gradual transition towards an elliptical 
cross section, with the 'effective' aspect ratio at ■ x *= 40 (y MflY /z MflY ) 
reduced from 4/1 to 2.3/1. The bulges due to wave/shear layer interactions 
are quite pronounced at x ■ 10 and 20. Note that the 4 contours at x * 10 
and 20 nearly coincide since the jet has fully expanded at x = 10 and then 
recompresses in going from 10 to 20. The predicted pressure variation along 
the jet axis of symmetry is exhibited in Figure 45. 
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FIGURE 32. U Contours For 4/1 Balanced Pressure Mixing Problem In 
”™~ Symmetry Planes . 
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FIGURE 33. 6 Contours For 4/1 Balanced Pressure Mixing Problem In 

Symmetry Planes. 
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FIGURE 35. 


K Contours For 4/1 Balanced Pressure Mixing Problem In 
Symmetry Planes. 













| FIGURE 38. Overall Jet Structure In Symmetry Planes As Exhibited By i 

Contours For 4/1 Balanced Pressure Case. 
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FIGURE 39, i Contours In Cross-Flow Plane At x * 60, 80, And 100 For 4/1 
Balanced Pressure Case. 
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FIGURE 41. Centerline Variation Of T; 4/1 Balanced Pressure Case. 
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FIGURE 43. 


Pressure Contours In Symmetry Planes For 4/1 Underexpanded 
Case. 
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3.10 Multiple Ondereipanded Square Jet Problem 


A calculation analogous to that of the underexpanded square jet problem 
(Section 3.7) was performed for an infinite array of square underexpanded 
jets as exhibited in Figure 46. The jets are initially 2x2 and their axes 
of symmetry have a spacing of 4. The lower portion of Figure 46 exhibits a 
schematic of the predicted shear layer growth and shock pattern for 0 < x < 
40 in the symmetry plane z = 0. The calculation was performed utilizing 
symmetry conditions on the planes z = 0, y = 0 and y = 2. The predicted 
pressure contours and i contours in the plane z = 0 are exhibited in Figure 
47. The flow is characterized by jet-induced bow shocks traversing the flow 
in the spanwise direction as readily evidenced by the pressure contours of 
Figure 47. The species contours expand /contract in accordance with the 
passage of the shock waves. The individual mixing regions coincide at x ~ 
35 as indicated by the coalescence of the outer 6 contours. 6 contours at 
the axial stations x = 0, 10, ..., 60 (Figure 48) exhibit the growth his- 
tory of the individual jets and their interaction for x > 35, with pronoun- 
ced wave interactions producing the corner region bulges in the early growth 
stage (0 < x < 20) . The mixing characteristics in the symmetry planes z ■ 0 
(TOP VIEW) and y = 0 (SIDE VIEW) are exhibited by the velocity and turbulent 
kinetic energy profiles at x = 0, 10, ..., 60 provided in Figure 49. 





FIGURE 46 . Schematic Of Multiple Square Jet Problem And Shock/Shear 
Layer Pattern In z = 0 Symmetry Plane. 
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FIGURE 49. Axial Velocity And Turbulent Kinetic Energy Profiles For 
Multiple Jet Problem. 
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4. NUMERICAL PROCEDURES FOB ANALYZING SUBSONIC MIXING 
PROBLEMS AND NUMERICAL STUDIES 


4.1 Overview of Pressure-Split PNS Subsonic Marching Approach 

One of the most popular procedures for the analysis of 3D subsonic 
mixing problems has been the SIMPLE procedure of Patankar and Spalding**. 
Numerous codes have been developed based on this procedure for both free and 
ducted jet mixing problems. In the SIMPLE procedure, the flow solution is 
obtained by a stepwise parabolic marching integration combined with a 
quasi-elliptic cross-flow integration. To permit spatial marching, the 
pressure field is split with the streamwise component imposed (or arrived at 
from global continuity constraints for ducted flows) and the cross-flow 
variation arrived at from the coupled solution of the continuity and 
cross-flow momentum equations. 

The SIMPLE procedure utilizes an iterative cross-flow integration 
approach whereby: 

(1) the cross-flow momentum equations are independently solved to 
yield the cross-flow velocities, V and W; and, 

(2) a pressure correction equation arrived at from continuity with 
'partial' cross-flow momentum coupling is used to upgrade the 
pressure field. 

Steps (1) and (2) are iterated upon until a converged cross-flow solution is 
obtained. This iterative cross-flow procedure is not always stable and 
according to Carroll**, 'most of the instabilities are caused by the coupl- 
ing between the equations'. To reduced these instabilities, significant 
underrelaxation is required between the iterative sweeps. In addition, in 
complex regions such as the initial corner regions of rectangular jets, 
convergence problems occur which sometimes cannot be remedied by under- 
relaxation. 

In analyzing 2D subsonic mixing layers, a very different 'philosophy' 
of approach is often taken for unbounded problems, which more closely follows 
standard boundary layer precepts. Thus, in the 2D boundary layer cross-flow 
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integration procedure of Bradshaw and coworkers 40 , 41 the cross-flow (normal) 
velocity is obtained from the continuity equation while the cross-flow 
pressure variation is arrived at from the normal momentum equation. Here, 
the roles of the cross-flow equations are reversed. This approach has 
better convergence characteristics than the SIMPLE procedure, and, if the 
pressure variation is negligible, reduces to a standard boundary layer 
solution requiring no iterations. 


Recently, Dash and Sinha* have reformulated the 2D cross-flow equations 
to obtain an independent equation for the cross-flow velocity which requires 
no assumptions about the cross-flow pressure or density variations. This 
provides a noniterative solution to the cross-flow equations for unbounded 
flows in situations with large cross-flow variations, and thus, improves 
upon the iterative procedure of Bradshaw and coworkers. An extension of 
this noniterative cross-flow procedure to 3D free jet flows is discussed in 
this report. The numerical framework for the parabolic (streamwise) portion 
of the jet mixing solution parallels that of SCIP3D as will be summarized 
below. 


In the new cross-flow approach formulated: 


(1) continuity-based equations for the cross-flow velocity components, 
V and W, are arrived at via combining the continuity and cross 
flow momentum equations, and, differential equation-of-state 
relations, subjected to fixed values of parabolic variables, 0 
(streamwise velocity) and H (total enthalpy) arrived at apriori; 

(2) the V (or W) equation is coupled to the W (or V) and P dependent 
variables only through vertical (z direction) derivative terms - 
thus in the absence of local cross-flow, no coupling occurs and V 
(or W) can be solved independently as in the 2D approach of Dash 
and Sinha; 

(3) recognizing that V and W must be closely coupled in vortical 
(e.g., corner) regions, the V and W equations are solved 
concurrently in a coupled fashion with the cross-flow pressure 
field lagged; 

(4) the cross-flow pressure field is updated using a cross-flow 
momentum equation - if the pressure changes are negligible, the 
solution is converged and no iterations are required; if the 
pressure changes are substantial. Step (3) is repeated using the 
updated pressure field; and. 
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(5) the direction of the cross-flow integration sweeps are alternated 
at subsequent steps to eliminate biasing of the difference 
relations. 


The present coupled, explicit cross-flow solution procedure is restricted to 
free (unbounded) jet flowfields with quarter plane symmetry. For other 
subsonic jet mixing problems with different types of boundary conditions, 
the cross-flow solution procedure requires modification. 


4.2 Parabolic Marching Procedure 


The 'parabolic portion* of the pressure-split PNS marching integration 
is performed using a variant of the SCIP3D integration procedure with the 
following modifications : 


(1) the streamwise (axial) momentum equation is solved with a - 0 so 
that Eg = pU 2 (rather than P+pU 2 ) and the pressure gradient, 3P/3x, 
is treated as a prescribed forcing function (see eq. 1); 

(2) an upwind variant of the MacCormack predictor/corrector algorithm 
is employed ( see ref . 1 ) ; 

(3) the wave portion of the jet growth rule given by eq. (19) is 
eliminated (i.e., dz^/dx = Zg(3f/3z) jj^/f a ) ; 

(4) a simple subsonic decode procedure is utilized to obtain the non- 
conservation variables from the E conservation array after both 
the predictor and corrector steps, viz. : 

f = E f /pU 

where : 

E f = pUf 

and; 

(5) the pressure field is approximated by: 

ap 3P 1 3P 3P (36) 

35 = H +a y3^ +a z3f 

which is solved in an ex plicit predictor/ corrector fashion using 
central differences, with 3P/3x specified. 


At the completion of the predictor/corrector 'parabolic integration', 
'firm' values of the variables 0, H and i> are obtained, and 'approximate' 
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values of the variables V, W and P are available which do not satisfy the 
continuity equation. The values of V, W and P are then revised to satisfy 
the coupled continuity and cross-flow (V,W) momentum equations with U, H and 
t fixed , and, with the approximate values of V, W and P used in evaluating 
coefficients for the cross-flow integration to be discussed below. No 
iterations between the parabolic and cross-flow sweeps are performed; 
rather, the axial step size is reduced in regions with large gradients (such 
as the initial jet mixing region) to minimize errors incurred by this 
approximation . 

4.3 Cross-Flow Equations 


Expanding the mapped form of the continuity equation yields the form 
listed below (eq. 37). 


|£ + v|£ + w|£+ ( 

3? 3t\ 3 


P b , 


3V 


+ ( 


p 2i 

u 


. aw 


- £ I - 4? + 4-r + a_ |? 1 


D 


3? y d-q z 3 C 


(37) 


In eq. (37), b , b . a„ and a_ are mapping transformation terms and the 
y z y z 

derivatives of the streamwise (axial) velocity, D_, U .and U are known from 

the parabolic integration performed apriori. The density derivatives in the 

continuity equation are replaced by derivatives of pressure and cross-flow 

velocities, as well as known derivatives of streamwise velocity and total 

enthalpy. This is accomplished by differentiating the perfect gas state 

relation (eq. 38) which yields the density derivative p given by eq. (39) 

a 

where a represents the 5 , ti or 4 directions. 


H = -2— - + i (o* + V* + W*) 
7-1 P 2 


which yields: 


P " p + (r-l) tUU + VV + WW - H ] 

a a 1 a a* a a a a 


(38) 


(39) 


- 83 - 



An important consideration concerns the treatment of the pressure 
gradient term, P^. Here, equation (36) is utilized which brings in the 
pressure-split approximation. Substituting the p Q expression into equation 
(37) with given by equation (36) yields the Modified Continuity Equation 
given by equation 40. 

Modified Continuity Equation 



a 


B = - - [y+_L]u +[-j£ r -yV]U 

P pU yM 2 5 yM 2 » 


a ~ 

[ - yW ] U + I [ H + VH + WH ] 

yM* t U l ri t 


(41) 


The cross-flow momentum equations in the y and z coordinate directions 
in transformed £, ti and t coordinates are listed below: 




(42) 


( 43 ) 
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where and are the parabolized stress terms, and y, p, V and W are 
defined by: 


y = ( y-1 ) / y 
p = pD* 

V = [b 1 - a ] 
y U y J 


W 


tb , J - “s’ 


The Modified Continuity Equation (eq. 40) and these two momentum 
equations comprise three equations for the cross-flow variation of three 
variables, V, W and P, with the parabolic variables 0 and H known apriori. 
By manipulations analogous to those of reference 8, continuity-based 
equations for V and V can be obtained which would be independent of all 
other parameters in the absence of cross-flow. The procedure for obtaining 
the 'W equation* is described below. 


Starting with the Modified Continuity Equation, we replace the terms 
and with their equivalents obtained from the V (eq. 42) and W (eq. 43) 
momentum equations, respectively. The term is given by: 


~v ~v \ b v ~ ~ 

(YfrJV* - (y^) I - -J P - W - WV ] 
0 5 U pU pu T) t\ $ 


The P term is given by: 

<b J)P fc - <*jp> l - W - VW - WW ] 
z P t P *1 t 


(44) 


(45) 


Substituting the expressions given by eqs. (44) and (45) into the Modified 
Continuity Equation, and grouping like terms, yields the expression: 



( 46 ) 
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Following the analogous steps, the 'V equation' listed below is obtained 



4.4 Cross-Flow Integration Procedure 

Referring to Figure 3, the cross-flow plane, 5 = constant, is sub- 
divided into IMAX-1 equally spaced intervals, An, where n = (1-1) / (IMAX-1) , 
and, JMAX-1 equally spaced intervals, A$, where ? = (J-l)/ (JMAX-1) . For 
quarter plane symmetry, the boundary conditions are as follows: 


W = 0 

on 

5 = 0 

(J=l) 

V = 0 

on 

n = o 

(1=1) 

^3 

II 

*0 

W 

on 

r 1 

( R = 1 

(J=JMAX) 

(I=IMAX) 


Equations (40) and (41) are solved in the following manner on the plane £ ■ 
constant : 

(1) the turbulent stress terms and are treated as forcing 

functions terms; 

(2) the pressure gradient terms and are lagged ; 

(3) the W equation (40) is integrated sweeping upward from 5 = 0 to 1 
( i . e . , from J = 1 to JMAX ) ; 

(4) the V equation (41) is integrated sweeping spanwise from q = 0 to 
1 (i.e., from I = 1 to IMAX); 

(5) the W and V equations are solved concurrently in a coupled manner 
which entails either: 


C-* 
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(a) sweeping upward from J to J + 1 with a nested spanwise sweep 
from I = 1 to IMAX; or, 

(b) sweeping spanwise from I to I + 1 with a nested upward sweep 
from J = 1 to JMAX . 

(6) after the coupled W and V solutions are obtained, the cross-flow 
pressure field is upgraded: 

(a) sweeping downward from £ = 1 to 0 (i.e., from J = JMAX to J = 
1) using the W normal momentum equation if an upward (a type) 
sweep was used in Step (5); or, 

(b) sweeping spanwise from n = 1 to 0 (i.e., from I = IMAX to (I 

= 1) using the V normal momentum equation if a spanwise (b 

type) sweep was used in Step (5). 

Referring to Figure 50, the V and W equations are solved using the 
triangular-type difference elements exhibited. The W equation is integrated 
upward in the £ direction from known values of W (=0) on the plane of 
symmetry, 5=0. The V equation is integrated spanwise in the ii direction 
from known values of V (= 0) on the plane of symmetry, n = 0. The upward/ 

spanwise integrations in the 5/n directions are performed using a two point 

(trapezoidal) procedure (as in the 2D cross-flow analysis of reference 8). 
The other family of derivatives (tj for W equation, 5 for V equation) are 
evaluated using central differences. 

4.5 Pressure-Split PHS Subsonic Square Jet Calculation 

This square jet calculation (Mach .5 into quiescent stream) was 
performed for a nozzle having the exit plane boundaries y = +1 and z = +1. 
The calculations were performed utilizing quarter plane symmetry with the 
initial grid extending to y = 2 and z = 2 yielding a substantial buffer 
region of uniform exhaust flow about the jet. Figure 51 depicts predicted 4 
contours at x = 20 obtained utilizing both the ke and kW turbulence models 
with the full pressure-split methodology. The faster mixing predicted by 
the ke model coincides with past experience for low speed axisymmetric jets. 
The pressure split solutions yield jet contours which have distinct kinks in 
the outer mixing region. This is to be contrasted with the smooth contours 
obtained by suppressing the effects of convective velocities on the solution 


- 87 - 





FIGURE 50. V And W Grids And Diagonal Coupling Nomenclature For Subsonic 
Cross-Flow Solution. 
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(i.e., by enforcing both V and W = 0 and thus generating a 'diffusion only' 
type solution) . 

The kinks in the 6 contours appear to be associated with the 'abrupt' 
breakdown of the initial vortex sheet surrounding the jet with the subse- 
quent formation of a pair of counter-rotating streamwise vortices 
symmetrically spaced about the corner region, as schematized in Figure 51*. 
A vector plot of the cross flow velocity field corresponding to the PNS/kW 
solution at x = 20 is given in Figure 52. We note that the jet entrainment 
velocities on the symmetry planes are in accord with values for planar jet 
mixing solutions. The 'flattening' of the outer 6 profiles along the 45° 
plane (y = z) is attributed to the significantly enhanced entrainment 
produced by the 3D corner mixing. The 'kinks' in the outer 6 profiles 
correspond to the positions where the streamwise vortices 'pump fluid out ' 
of the jet. 

The subsequent development of the jet out to x = 100 (PNS/kW solution) 
is exhibited by the 6 = 0.1 contours at x = 40 (inner contour is 6 = 0.9) 
60, 80 and 100 superimposed on the cross-flow velocity vector field (Figure 
53). The corner-induced streamwise vortex pattern persists in the cross- 
flow plane despite the nearly circular behavior of the streamwise parabolic 
solution. The variation of i and turbulent kinetic energy, k, along the jet 
centerline are exhibited in Figure 54, and, are compared with 'correspond- 
ing' planar and axisymmetric jet solutions (viz., solutions having nozzle 
exit heights/ diameters of 2) predicted using the SPLITP jet mixing model. 

4.6 Pressure-Split PHS 4/1 Aspect Ratio 

A 4/1 aspect ratio jet was analyzed having conditions identical to 
those of the square jet case. The nozzle exit extended to y = +2 and z = 
jK).5, while the initial grid, employing quarter plane symmetry, extended to 
y = 2.5 and z = 1 (i.e., a uniform buffer of freestream flow of width ~ 0.5 


The numerical approach may also have some influence on the vortical 
behavior predicted, e.g., the spatial marching procedure provides no 
mechanism for streamwise pressure gradients and/or upstream influence 
effects to be accounted for, which might influence the results. 
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FIGURE 53. Subsonic Square Jet Calculation; t = 0.1 Contours Super 
imposed On Cross-Flow Vectors Maps At x = 40 - 100. 







surrounded the jet initially) . Figure 55 depicts predicted jet contours at 
x = 20 using the kW turbulence model with the pressure split, PNS approach, 
and, with the 'diffusive limit' approach whereby the convective velocities 
are suppressed (V = W = 0) . As in the square jet case, the PNS 6 contours 
exhibit pronounced kinks associated with the corner generated streamwise 
vortices depicted in Figure 55, while the case with convective velocities 
suppressed yields smooth, kink -free contours. 

The development of the jet out to x = 100 is exhibited by the 6 = 0.1 
contours in the cross-flow plane at x = 40 (inner contour is i = 0.9) 60, 80 
and 100 (Figure 56) with a schematic of the streamwise vortex pattern 
superimposed. By x = 100, the jet contour takes on an elliptic appearance 
with the vortex-induced kinks no longer evident. The 4> contours in the y = 
0 and z = 0 symmetry planes out to x = 300 are depicted in Figure 57, 
indicating that by x = 300, the flow is nearly circular. Profiles of i and 
turbulent kinetic energy, k, in the two symmetry planes at x = 100, 200 and 
300 are exhibited in Figure 58. The variation of axial velocity, U, and 
turbulent kinetic energy, k, along the jet axis of symmetry are depicted in 
Figure 59. 

4.7 Parabolic Marching Approach (P = P(x) Only) 

An alternate procedure for analyzing subsonic jet mixing problems 
involves the solution of a reduced set of parabolic 'boundary layer 
equations ' . In this approach, the cross-flow grid must be aligned so that 
one coordinate direction is essentially normal to the jet boundary. This 
can be accomplished using the hybrid adaptive rectangular/ cylindrical grid 
discussed previous (see Figure 4). Utilizing this adaptive grid, the z 
(cartesian) and r (cylindrical) directions align themselves normal to the 
jet boundary. Thus, W represents the velocity in the normal (entrainment) 
direction and V represents the cross-flow velocity. The equations solved 
neglect pressure variations in the cross-flow plane (P = P(x) only). The 
solution procedure employs a single parabolic solution sweep whereby the V 
component of velocity is determined from the V-momentum equation while the W 
component of velocity is determined from the continuity equation 
aposteriori. The W-momentum equation is neglected consistent with the 
approximation of constant pressure in the cross-flow plane. 
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X = 40 


X = 80 



FIGURE 56. Rectangular Subsonic Jet Calculation, Contours At x = 40, 60, 
80 And 100. 
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4.8 Parabolic Rectangular Jet Analysis in Adaptive Cylindrical Coordinates 


The calculation described here utilized an adaptive cylindrical grid to 
analyze a 2/1 Aspect Ratio rectangular jet (M^ = 1, Tj = 800 °K) exhausting 
into a Mach .3 external stream (Tg = 278 °K) . The cylindrical grid employed 
21 points in the radial direction and 11 4 planes (0 < 6 < n/2) . The cal- 
culation was performed using the simplified parabolic boundary layer run 
option described above. Figures 60 and 61 exhibit predicted 4 and T 
contours in the y = 0 and z = 0 symmetry planes. Figures 62 and 63 exhibit 
predicted 4 and T cross-flow contours at x = 5, 10, 20 and 30. 4 and T 
cross-flow contours at x = 50 (Figure 64) show a nearly circular behavior. 

4.9 Comparison of Parabolic and Supersonic PHS Predictions for 3/1 Aspect 

Ratio Rectangular Jet Osing Hybrid Adaptive Grid 

This calculation was performed for a hot Mach 3 exhaust stream (T^ = 
2000 °K) exhausting into a colder Mach 3 external stream (Tg = 500 °K) to 
compare the simplified parabolic run option with a full supersonic PNS 
calculation. The calculations were performed using the rectangular/ 
cylindrical hybrid adaptive grid with a moving origin. Figure 65 compares 
cross-flow 4 contours for the two solutions at x =20. Both exhibit the 
same 'kinked' behavior in the corner region. Figure 66 depicts the cross- 
flow velocity vector field at x = 20 for the parabolic solution - note that 
the cylindrical origin (here at y = 1) moves inward to maintain the same 
dimensions on the 4 = 0 and 4 - n/2 planes. Figure 67 exhibits the 
evolution of the jet boundary contours from rectangular to circular. 
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FIGURE 60 


i Contours For 2/1 Rectangular Jet, Parabolic Run Option 
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FIGURE 64. t And T Cross-Flow Contours At x * 50 For 2/1 Rectangular Jet 
Problem. 
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SUPERSONIC PNS SOLUTION 



FIGDBE 65, 6 Cross-Flow Contours For 3/1 Rectangular Jet Problem, x = 20 

Supersonic PNS vs. Parabolic Solutions. 
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CROSS-FLOW VECTORS 



FIGURE 66. 


Cross-Flow Velocity Vector Field At x = 20 For 3/1 Rectan- 
gular Jet Problem, x = 20 Parabolic Run Option. 
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5. SPHERICAL PROCEDOBE FOB ANALYZING SUPERSONIC JETS EXHAUSTING 
INTO SUBSONIC EXTERNAL STREAMS AND NUMERICAL STUDIES 


5.1 Review of 2D Numerics for Supersonic/Subsonic Jet Mixing Problem 

The spatial marching PNS analysis of a supersonic jet exhausting into a 
subsonic external stream requires: 

(1) switching the numerical algorithm from full PNS at supersonic jet 
grid points to parabolic or pressure-split PNS at subsonic jet 
grid points; 

(2) matching the supersonic and subsonic PNS solutions at the jet 
sonic line; and, 

(3) providing for a global pressure iterative methodology to account 
for the overall interaction of the jet and subsonic external flow. 

At the 2D level, the above procedures have all been made operational. The 
analysis of the mixed supersonic/subsonic problem in the absence of global 
external flow interactions (e.g., the external pressure field is assumed to 
be constant) has been performed with the SCIPVIS model 1 , and, numerous 
comparisons with the data of Seiner have been performed which demonstrate 
the viability of this approach (see refs. 2-5). A typical solution is 
exhibited in Figure 68 (from ref. 1) which shows the ability of SCIPVIS to 
accurately predict the complex jet pressure field. A key feature of SCIPVIS 
was the characteristic-based methodology developed to couple the supersonic 
and subsonic solutions at the jet sonic line as schematized in Figure 69. 
The coupling (at the grid point I*) is performed with the outer jet boundary 
pressure, P g (x), prescribed and a viscous-characteristic constraint just 
under the jet sonic line (along X+) relating the supersonic pressure and 
flow angle. In SCIPVIS, the normal pressure variation across the subsonic 
outer layer is neglected and P e (x) is imposed at the sonic line yielding the 
flow angle there. Then, the variation in flow angle across the jet is 
obtained from the continuity equation yielding the entrainment (suction) 
velocity, V e (x), at the jet outer edge. 
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JET MIXING LATER 



FIGURE 69 


2D Supersonic/Subsonic Coupling At Jet Sonic Line 





5.2 gonadaptlye Rectangular Matching Boundary Approach 


For 3D jets, the coupling problem is complicated by the shape taken by 
the sonic line in the cross-flow plane, as well as by initial corner region 
vortical effects. In our initial work, mapped rectangular coordinates were 
employed and sonic line coupling was performed along the 'largest' embedded 
rectangle that could fit within the sonic line in the cross-flow plane (Fig. 
70) . The pressure along the matching boundary was lagged yielding V and W 
at all I* and J* matching points which served as inner boundary conditions 
for the subsonic coupled V-W cross-flow solution. The subsonic cross-flow 
solution was swept outward from the matching boundary to the outer boundary, 
with y and z sweeps independent of each other until the I*, J* matching 
point. The nested subsonic corner region (I > I*, J > J*) was treated using 
specialized coupled logic. 

This approach was applied to the analysis of an underexpanded square 
jet exhausting into a subsonic stream. Reasonable results were obtained for 
the initial jet expansion (since the square grid coincided with the initial 
sonic line shape) , but became quite poor when the reflected waves (which 
recompress the jet turning the stream surfaces inward) reached the now 
quasi-circular sonic line. The square representation of the sonic line 
artificially delayed the recompression waves from interacting with the 
corner region sonic line as shown in the sequence of cross-flow contours 
displayed in Figure 71. As exhibited, the initial expansion of the jet is 
well-predicted, but the recompression process (x > 4) leads to highly 
unrealistic cross-flow contours. 
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3.3 Variable Rectangular Grid Matching Boundary Approach 


To remedy this problem, a variable matching boundary which adapted 
itself to the sonic line cross-flow variation was implemented, still working 
in mapped rectangular coordinates (Figure 72) . Using this approach, the 
matching logic at the sonic line is quite complex, and the solution of the 
coupled V-W subsonic cross-flow equations entails very cumbersome logic. 
Results using this adaptive approach are exhibited in Figure 73. The solu- 
tion in the corner region ultimately goes unstable due to complexities in 
concurrently handling the initial subsonic corner-induced vortices and the 
reflected compression wave system. 

3.4 Fully Adaptive Matching Boundary Approach 

The current approach for dealing with this complex problem involves 
utilizing adaptive coordinates with the cross-flow directions aligned: 

(1) along the jet entrainment (inflow/outflow) direction, essentially 
normal to the sonic line; and, 

(2) normal to this direction and thus aligned with the sonic line and 
also normal to the plane/axis of symmetry. 

In this approach, a parabolic procedure is employed for analyzing the 
subsonic region of the flow which neglects the momentum equation in the 
entrainment aligned direction. Results for an 8/1 velocity ratio balanced 
pressure square jet problem (Mj = 2, Tj « 2000 °K; Mg = 0.5, Tg = 500 °K) 
performed with adaptive cylindrical coordinates are exhibited in the next 
several figures. The calculations were performed using a fully parabolic 
approach and the adaptive supersonic/ subsonic PNS approach discussed here. 
Parabolic vs. PNS solution contours of 4 and T in the jet symmetry plane are 
exhibited in Figures 74 and 75, respectively. Cross-flow contours of 4 and 
T at x = 25 are compared in Figures 76 and 77. The 'crude' computational 
grid employed to perform this calculation is exhibited at x = 25 in Figure 
78. Pressure contours at x = 25 are exhibited in Figure 79. Note the close 
similarity in 4 contours showing that mixing processes are nearly identical, 
but, variations in T contours indicating the influence of the variable 
pressure field in the PNS case. 
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MATCHING POINTS 
VARIABLE 1* AND J* 


FIGOBE 72. Variable Rectangular Grid Matching Boundary. 
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t CONTOURS 


CROSS VELOCITY VECTORS 



FIGURE 73 


Species (6) Contours And Cross-Flow Velocity Vectors For 
Underexpanded Square Jet Into Subsonic External Stream Using 
Version II Supersonic /Subsonic Coupling (Variable Grid 
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FIGURE 76. 6 Contours At x = 25 For Balanced Pressure 

Supersonic Square Jet Problem. 
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FIGURE 78. Adaptive Cylindrical Grid In Cross-Flow Plane At x = 25 For 
Balanced Pressure Square Jet Problem. 
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FIGURE 79. 


Perturbation Pressure Contours In Cross-Flow Plane At x = 25 
For Balanced Pressure, Square Jet Problem. 



Preliminary results for an underexpanded subsonic/supersonic square jet 
problem with conditions comparable to those of the Seiner case 1 -* (see Fig. 
68) are exhibited in Figures 80 and 81. Figure 80 exhibits i and P contours 
at x = 4. The predicted variation of pressure along the jet axis for the 
first shock cell is exhibited in Figure 81. To check out the code's 
performance, an attempt was made to analyze the above problem in the 
axisymmetric limit and compare it directly with SCIPVIS results. The 3D 
code utilized the same grid as SCIPVIS in each cross-flow plane, i = 
constant. The 3D calculation was performed utilizing quarter plane symmetry 
employing five i planes. Figure 82 compares SCIP3D predicted shear layer 
boundary contours (tf = 0.05 and 0.95) in the symmetry planes with corre- 
sponding SCIPVIS results. Figure 83 compares predicted centerline pressure 
variations. The comparisons indicate that further work is required to deal 
with near-axis effects in the 3D code to get the two solutions to agree more 
closely. 



0.0 0:1 0.4 0.6 0.8 1.0 
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FIGURE 80. Contours Of 6 And P For Underexpanded Subsonic/Supersonic 
Square Jet Problem. 
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FIGURE 82 


Comparison Of SCIP3D And SCIPVIS Predictions Of Shear Layer 
Boundaries For Underexpanded Axisymmetric Subsonic/Supersonic 
Jet Problem. 




Pressure ( P/P ( 
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FIGURE 83. 


Comparison Of SCIP3D And SCIPVIS Predictions Of Jet Center- 
line Pressure Variations For Underexpanded Axisymmetric 
Subsonic/Supersonic Jet Problem. 



6. SUMMARY 


The SCIP3D 3D PNS jet mixing model and its application to a variety of 
rectangular jet mixing problems have been described in this report. For 
fully supersonic free jet problems (supersonic jet into supersonic external 
stream) , the model has been demonstrated to predict plume flowfields for 
mildly underexpanded cases, and can be used to support plume related 
engineering studies at a research code level. The SCIP3D model also has 
current applicability to supersonic ducted problems (e.g., 3D nozzles/ 
diffusers and ducted jet mixing problems) but lacks the ability to resolve 
the near-wall boundary layer due to the explicit numerics utilized. 

An extended version of SCIP3D (entitled SCHN0Z3D) has been developed 
for applications to scramjet propulsion problems. SCHN0Z3D contains 
generalized finite-rate chemistry capabilities and has been applied to 
scramjet combustor and nozzle/ plume flowfields. Features and applications 
are to be reported in the Fourth National Aerospace Plan Tech Mat Symposium, 
Monterey, CA, February 1988 and the 17th JANNAF Plume Technology Symposium, 
NASA/LRC, Hampton, VA, April 1988. Exploratory work to extend this model to 
analyze the near-wall boundary layer is in progress utilizing a time- 
iterative procedure. 

With regard to the analysis of supersonic jets exhausting into a 
subsonic external stream, the current noniterative methodology requires 
gridding which adapts to the sonic line shape in the cross-flow plane. 
Successful runs have only been performed by invoking simplifying assumptions 
for the cross-flow pressure variation in subsonic regions (e.g., P = P(x) 
only) . It appears that the noniterative methodology is quite cumbersome to 
deal with and that efforts should focus on time-iterative methodology on a 
plane to plane basis (e.g., in a spatial marching mode of operation with 
upstream influence effects suppressed) . Other recommended upgrades to 
SCIP3D would include casting the equations in generalized coordinates and 
incorporating adaptive procedures to pack grid points into shear layer 
regions . Upgrading turbulence models from the two-equation level to the 
multi-scale or algebraic Reynolds stress level would also appear requisite 
to deal with vortical aspects of 3D mixing. This should be followed by a 
detailed validation study using the available rectangular jet mixing data 
base. 
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